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ABSTRACT 


The widespread availability of inexpensive high-speed computers has led to 
the development of complex, detailed, technical models of combat. These high 
resolution computer simulations and wargames are touted by their proponents as low- 
cost alternatives to extensive, high-cost field training exercises for the training of 
combat leaders. The validity of these simulations as models of combat, and thus as 
useful training tools is unproven. Direct comparison of simulations with field training 
exercises 1s often frustrated by the inherent complexities in each, and the shortage 
of quality data from field exercises. This thesis examines the feasibility of comparing 
these systems indirectly through the use of surrogate analytical models. A discrete 
mathematical surrogate model is examined, and used to generate probability surfaces 
for comparison. We consider several techniques of fitting surrogate models to 
combat data, including simulated annealing, steepest descent, and multiple regression. 


Areas for further research are discussed. 
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THESIS DISCLAIMER 


The Reader is cautioned that computer programs developed in this research may not have 
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not be considered validated. Any application of these programs without additional 
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I. INTRODUCTION 


Comparing two high resolution models of combat has proven to be a difficult 
task. This thesis uses a simple discrete dynamical model of combat as a surrogate 
for the complex models being studied. Our surrogate model is fit to simulated battle 
data generated by a mathematical model, and used to generate probability surfaces. 
These surfaces are compared using a randomization test [Ref. 1] which proves 
effective at identifying battle data sets from similar simulations models and at 
distinguishing between battle data sets from different simulation models. 

Land combat has evolved into a Series of complex, combined arms battles fought 
with extremely expensive high-technology weapons. Many recent efforts to model 
this process have relied on computer simulations. These simulations attempt to 
capture the processes of combat by simulating the detailed interactions between 
individual combat elements (e.g. tanks, armored personnel carriers, and artillery 
pieces). As a result, these high resolution simulations are highly complex and 
expensive in their own right. 

Having paid the price to develop these simulations, users rightly want to know 
if they replicate the combat processes modeled. This question has to a large extent 
gone unanswered. A great deal of data analysis, experimentation, and field testing 
goes into developing procedures used to simulate the various "microscopic" combat 
processes, such as searching for targets, identifying targets as hostile, and engaging 
hostile targets. Verifying that these procedures adequately model the processes for 
which they are designed is relatively straight forward. The interactions amongst all 
the various underlying processes are not well understood, and that makes the 
validation of the overall model difficult. 

To validate a high resolution combat simulation such as Janus(A) [Ref. 2], an 
army battalion task force air-land combat simulation, one might try comparing some 


measure of effectiveness (MOE) from an actual combat engagement to the same 


measure of effectiveness for a simulation of that battle in Janus(A). However, not 
much data are available on actual engagements of high-tech armored forces, and 
actual combat data tends to be clouded by haphazard data reporting and collection 
techniques. The next best procedure might be comparison to data collected from 
realistic field training exercises such as those conducted at the United States Army’s 
National Training Center (NTC) at Fort Irwin, California. The data from NTC, 
although collected more systematically than combat data, also has its shortcomings 
[Ref. 3], and it is expensive to collect. For these reasons it is difficult to get the 
desired number of replications of similar battles from NTC. 

How then can we validate a model such as Janus(A) if there is no good source 
of data for comparison? This thesis attacks that problem by exploring ways to use 
analytical (mathematical) surrogate models of high resolution simulations, such as 
Janus(A) and field exercises at NTC, fitted to small battle data sets. Chapter II 
outlines the methodology proposed to compare two combat models. Chapter III 
examines the utility of simulated annealing, steepest descent, non-linear optimization, 
and multiple regression techniques for fitting the parameters of an analytical 
surrogate model to the simulation model’s data. Chapter IV presents an example of 
the methodology to compare similar and dissimilar battle data sets generated from 
a mathematical model, and finally Chapter V concludes this study with a discussion 


of the utility of the methodology and proposals for further research. 


Il. METHODOLOGY 


A. INTRODUCTION 

This chapter examines the methodology developed for comparing two high 
resolution models of combat (referred to as simulation models throughout this 
thesis). First, we examine the assumptions necessary for our method, and the 
Surrogate models considered (including a development of a discrete Lanchester 
model). We then discuss the techniques investigated for fitting battle data to the 
surrogate model and the measures considered for comparison. Finally, we examine 
the actual procedure used for comparing the two simulation models (i.e. two machine 
generated battle data sets representing simulation models). 

Throughout this chapter we will refer to two critical items which are defined 
below. 


e Battle trajectory. A battle trajectory is a numerical record of battle expressed 
as the number of each weapon system still active in the battle at the end of 
each time increment. 


e Battle data set. A battle data set is a collection of battle trajectories from the 
Same simulation model. 


B. ASSUMPTIONS 
Certain assumptions are necessary to implement this methodology. 

1. Battle trajectory data are obtained as the raw number of each weapon 
system actively involved in the battle, recorded at the end’ of each time interval, as 
shown in Table J; all time intervals are of the equal length. 

2. It is assumed that each weapon system involved in the battle is capable of 
engaging any other weapon system involved in the battle. This assumption is 
required when using Lanchester-type equations as the surrogate model (Lanchester 


Models of combat are discussed in Section D). 


TABLE 1. SAMPLE BATTLE TRAJECTORY. 


Battle Trajectory from NTC 
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3. The battle trajectories are assumed to be non-increasing (i.e. no 
reinforcements are allowed). This is often the case in computer simulation models. 
This assumption forces the combat processes’ mathematical representation to be 
convex, thus simplifying the fitting of the analytical model. 

4. Individual battle trajectories are assumed to be independent of all other 
battle trajectories. This is not a problem with computer simulations that are 
restarted using different random number seeds. However, this can be a problem with 
data obtained from field exercises such as those conducted at the NTC. Similar 
battles conducted by the same unit (or, even with the same opposing force) can 
display learning curve trends which could result in correlation between different 


battle trajectories. 


C. THE SURROGATE MODEL 
Choosing the surrogate model is a critical step in the analysis process. The 


model chosen must represent the underlying processes of interest. Ideally, the model 


should also be simple enough (i.e. have a small number of parameters to be 
estimated) to allow good fits to small data sets. 

The model chosen for our initial trials was a discrete analog of the Lanchester 
aimed fire model. This model was chosen because it is an attrition model; because 
it has a small number of parameters to be estimated; and because the weapon 
systems considered were all direct fire weapons (i.e. tanks, armored personnel 
carriers, and anti-tank weapons). A short introduction to Lanchester-type models is 


provided in the next section. 


D. LANCHESTER MODELS OF COMBAT 

In the late eighteenth century Antoine-Henri Jomini wrote extensively on what 
he called the principles of war [Ref. 4]. One of these principles was that an army 
should mass (concentrate) its forces when attacking. Over a century later, in 1914, 
a British engineer and inventor F. W. Lanchester set out to prove a hypothesis 
similar to Jomini’s [Ref. 5]. Lanchester conjectured that in "modern warfare" the 
concentration of forces was an effective tactic. To prove his hypothesis Lanchester 
developed a mathematical model of "modern warfare." Lanchester argued that with 
modern weapons it was possible for any combat element to engage any other on the 
battlefield, thus the attrition of a force should be proportional to the size of the 


opposing force. Lanchester presented the following model: 


rae = —q > a>o, x0) = Xo: (1) 
oo} aa . = 
me = myo; “D>O, Y(O) = Yo. (2) 


where 
X(t) = X force’s strength at time t, 
Y(t) = Y force’s strength at time t [Ref. 6]. 


This model can be solved by dividing equation (1) by equation 


(2): 





ax 
dt _ -ayY 
“ay px! ey 
ae 
Cx aad 
Tere | Seok (4) 
-bXdX = -aYdy. (5S) 
We then integrate both sides: 
E 
[ -bxdx = f[ -a Ydy (6) 
0 0 
E if 
= 12) 2 = -2y?2/, (7) 
Z 0 2 0 


-b(X*(t) - X?7(0)) = -a(Y¥?(t) - ¥2(0)). (8) 
Thus we have: 
Di Ke" XX" (eC) saa ee een, eu 


which is known as Lanchester’s Square Law. 

Lanchester’s model is a continuous approximation of the discrete combat attrition 
process. This formulation presents two problems for direct use as a surrogate model. 
First, in order to use Lanchester’s model to simulate combat on a computer, it is 
necessary to approximate the differential equations with difference equations. 


Second, the observations we are using to estimate the parameters in our model are 


taken at discrete time intervals. It thus makes more intuitive sense to formulate the 
model with difference equations from the start. Lanchester’s model expressed as a 


system of difference equations is: 


Xn — X, = -A@Y,3 a0, Y= YY, (10) 

ene be0, XK. = xX, (11) 
By iterating equation (10) one step we have: 

epee = ee Ye, (12) 


We then subtract equation (10) from equation (12) which yields: 


ane a Ker emma (Vim Yo), (13) 
and by substitution we get, 
Kim» = Bees Oe Oe (Slope (14) 
Thus, 
hae i= 2 Gea ual: (15) 


which is a second order discrete dynamical system. It can be solved using standard 
techniques for second order difference equations [Ref. 7]. Its characteristic function 
is given by: r* - 2r + (1 - ab) = 0. 

This equation has two real roots, 


r= ri>2) +4 -4(1 - aby 


2 


— Leave eye (16) 


Thus the general form of the solution to the discrete dynamical system for the X 


force is given by: 


X, = C,(1 + Jab)* + C,(1 - fab)”, (ie7e) 


n 


where C, and C, are constants determined by the initial conditions. Applying our 
initial conditions: 

Xy=% and X,-Xy = -ayo 
we obtain the following two equations in two unknowns: 


X= Xo C,(1 + fab)® + C,(1 - fab), (18) 
X, =X) - AVg = C(1 + ab)’ + C,(1 - Vad)’. (19) 


Solving for Cy and C, we get: 


6 aes 

C= = : te Yo - 
Thus the particular form of the solution for the X force is: 
= 7 = Yo = ef = Yo 


and similarly for the Y force: 


(1 = Jab)” , 425 














ne (1 + Jab)” 
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Yo b 
fe) 
D ey 


(1 = fab)”. (2 
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The graph of these solutions with two different values of Vab is shown in Figure 1. 


Y, = (ay ab) aes 


n 





From this graph we see that the solution trajectories are approximately linear over 


the domain of interest (the first quadrant) with slopes determined by the value of 
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Figure 1. Solution Trajectories. 


Vab. The larger the value of Vab the more intense the attrition process and thus the 
steeper the slope. We can interpret the Vab as a measure of the pace or intensity 
of battle (see Section F). 

Another value of interest is the time at which the losing force would be 
annihilated. We are unable to calculate this value for the continuous Lanchester 
model because the solution trajectories go to zero asymptotically. 

For the discrete model the annihilation time of the losing force is calculated by 
setting X(n,) = 0 and solving for n,. Thus we have: 


X(n,) =0 [2 - ; <M 


(1 eye) et 























Now provided the right hand side of this equation is positive, which we show below, 


we can take logarithms on both sides: 
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Similarly from equation (21) we get the annihilation time for the Y force: 


We now show that the quantities in the above equations are positive. We know that 
if Vab > 1 then at least one of the attrition coefficients a,b > 1. This implies that in 
the battle, at least one of the forces would be annihilated in one time period. Thus 
we are only interested in the case when Vab < 1, and we have adjusted the length 
of our time periods to insure that a,b < 1. Thus in order to show that we can always 


determine the battle termination time n, we need only show the following. 


Lemma. 
Assume X,Y, > 0 and 0<a,b< 1. 
x b Yo 
[ = - 2 <0 then — x, -— >9Q. 
y {2 0 ~ "9 {2 a 
Proof: 
a Xo 
Assume = ee () 
ie %0 5 
x 
Then se an 
p>” 9 
SO 
b Yo 
pie, > eae" >—, 
. 0 Yo 
ea ee 0. 
a 2 
Q.E.D. 


a, 


Using the formulas we have derived for the battle termination time we can now 


develop conditions for determining which force will be victorious. Settingn, = hy 


we get: 
Dies Xo b Yo 
Boge call Pspanotlo 
in b 2 In a Zz 
Xo a Yo {2 
—+,|-—y — +, |—x 
Sea? Bal ie a °| 
G Xo b Yo 
— _->_ — — xX. - — 
oe es 2 
Xo a Yo b 
— = — x 
Bs |e, ae 
SO 


a Xo 1] Yo b Xo a b Yo 

—y--li- —x,| = }|J—+,A- —X%,- —|}. 
{> s|(3 {2 3 KE ; 
Multiplying out the terms 


2 2 
a Yo XyYo | b Xo b 2 = *o% G2 
—— +X -~— |—-—_ = | —x,* - — +x - j= ; 
io Mare a (2s 4 00 [2x 
a b 
{230 ‘os one 


SO 
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fiiatels: 


Therefore if cA then the battie is a stalemate, 


= 
oi 
& | 


else if —— > 


then the Y force wins, 


else if An then the X force wins. 


This analysis of our model suggests a number of natural measures of battle which 


we will consider in Section F. 
1. Non-Homogeneous Lanchester 


Our model up to now has considered only homogeneous forces, such as a 
battle of tanks versus tanks. The models we wish to compare contain not only tanks, 
but armored personnel carriers and anti-tank weapons (they may also contain indirect 
fire weapons such a artillery; however, we will not consider them here). To include 
the differences in the attrition potential of the various weapon systems, we assign 


each weapon a weight, which we call its "fire power index." The total attrition 
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potential of the X and Y forces is calculated as a linear combination of the individual 


fire power indices [Ref. 8]. 


Vey J. 
i 


X= Phx, 02) 
j 


Where; 
Y = the total fire power of the Y force, 
y, = the number of Y systems of type 1, 


and 


f,, = the fire power index of a Y system of type i, 


and similarity for the X force. 

2. Stochastic Lanchester 
Our model has been deterministic thus far, with the battle results being fixed 
by the initial force levels on each side and the fixed attrition coefficients. There is 
intuitively much more to who wins a battle than simply the initial force levels. 
Leadership, training, momentum, and terrain are just a few of the critical factors in 
determining the outcome of any real battle. Introduction of these factors explicitly 
into our model, even if feasible, would make it too complex to fit to small data sets. 
We take a simpler approach by adding these factors as stochastic noise. Gomg back 
to the difference equation models, let us examine how this could be done. One 


simple approach is to add a random term to each equation: 


Ani] Ae Og te 
Ve tn = (OX eee (23) 
where Z, and Z, are random variables with given distributions whose parameters are 


to be estimated. 
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This approach, however, suggests that the noise is independent of the force 
sizes involved. This is intuitively unappealing as history shows the more complex the 
battle the more confusion likely to be present. We can incorporate this idea into our 
model by postulating that some factors, such as combat readiness and momentum, 
may be dependent on force size; while other factors, such as leadership and training 
may be independent of force size. Thus we separate these factors into two random 


terms: 


Xnel _ Xx, 
nae as Yr, 


-A(n) Y, + Z(n), 
- B(n) X, + Z,(") : (24) 


where A(n) and B(n) are random variables representing the force size dependent 
factors, and Z,(n) and Z,(n) are random variables representing the force size 
independent factors. 

Further we assume that these factors are predetermined at some level at the 
beginning of each battle based on the training and maintenance preparation of the 
forces involved. These values then vary about these fixed values based on noise 
introduced by the confusion and stress of combat. We therefore postulate that for 
any given battle the random variables (A(n), Z,(n), B(n), Z,(n)) can be represented 
by a constant plus some error (noise) random variable which is independent, 
identically distributed normal with mean zero and variance o%. 

Thus we write: 
-A(n) = -a + 6e,(M), 
Z(m) = 2, + €2(M), 
-B(n) = -b + €,(n), 
Ze = Z, + €; (n). 5) 
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Substituting these expressions into equation 24 yields the following system of 
equations: 
Xiu —%, =[-a + e,@]Y, + (z, + e-@I, 


Ve eg > Le, oe Sous €z (”)]. (26) 


n+1 


It should be noted that in these equations the a, z,, b, and Zy are only constant over 
the duration of battle. They will vary from battle to battle based on the combat 
preparation, force mix (e.g. the ratio of tanks to APCs), and the commander’s 


concept of operations. Thus over the battle data set we will have a set of values 
{4, 2,9 5, zy} for each battle which can be considered to be realizations of the 
random variables A, Z,, B, and Z,, where A, Z,, B, and Z, are jointly distributed with 
parameters to be estimated. 

Finally, we include the individual attrition potential of each weapon system 
(fire power indices) from equation 22. 
This completely specifies our model as follows: 

Xu. - X, = [-a + el, + bz, + e, 1, 


Yiu - ¥, = [-b + eX, +k, +e, 27 


where, 


x; = dA Xi rag 2 SM 


J 
Si, = the fire power index of X force weapon system type J, 
hy = the fire power index of Y force weapon system type j, 


X,, = the number of X weapon systems of type J active in the 
battle at the end of time increment 1, 


16 


Y, = the number of Y weapon systems of type j active in the 
battle at the end of time increment 1, 


e,{n) is iid Normal(0,o,”), 

é7 (n) is iid Normal(0, Oo,” Mr 

€,(n) is iid Normal(0,o,’), 

€z (n) is iid Normal(0,0,”), 

(note €,(%), €7 (7), €,(n), and €, (x) are also independent of each other) 


4, Z 5, 2, | is a realization of [A, Z,, B, Z,) 


and a,b > 0. 


E. FITTING TECHNIQUES 

A number of techniques are evaluated in this thesis for estimating the parameters 
of a system of differential/difference equations from small battle trajectories. Using 
the non-linear programming solver MINOS 5.2 with the General Algebraic Modeling 
System (GAMS) to do a simultaneous least squares fit on the data was found to be 
the most effective. However, suppose one is willing to make the simplifying 
assumption that in an individual battle the coefficient random variables in one 
equation are independent of the coefficient random variable in the other equation 
(i.e that [A, Z,] and [B, Z,] are independent, but A and Z, are still considered to be 
dependent, as are B and Z,). In this case the parameters can be estimated by 
standard regression techniques. The validity of this assumption can be tested by 
fitting a small number of battle trajectories, using the GAMS model which considers 
each equation simultaneously, and comparing the resulting estimated coefficients with 
estimates obtained from the standard linear regression fits of the same battle 


trajectories (i.e. those obtained by considering each equation separately). If the 
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comparison is favorable (i.e. the difference between the parameters fit by the two 
models is small), the savings in the computation time required for parameter 


estimation is great. 


F. MEASURES OF BATTLE 
In this section we consider measures derived from the difference equation model 
with only force size dependent terms; that is, all factors are considered to be 


dependent on force size. With this assumption equation 27 becomes: 


XxX. - %, = l-a4+ eI] ¥, . 


Y.,-Y,=[-b + e,(n)]X, . (28) 


n 


In defining the measures of battle the following definitions will be used. 


AX, = the change in the fire power of the X force during 
the time increment [n, n+1]; so AX, = X,,, - X,. 


AY , = the change in the fire power of the Y force during 
the time increment [n, n+1]; so AY, = Y,4,- Y,- 

X, = fire power of the X force at the beginning of the 

time increment [n, n+1]. 


n 


Y, = fire power of the Y force at the beginning of the 
time increment [n, n+ 1]. 


N =the number of time increments in the battle 


trajectory. 
N-1 A X. 
G = estimate of a for a single battle; io «YY, 
av=te 
N 
oa) 
b = estimate of b fora single battle; | on oe 
pa 
N 
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G, = the estimated value of a on [n, n+1]; 4, = : 
n 

: , SNE 

b_ = the estimated value of b on [n, n+1]; 5, = = 


1. Pace of Battle 


The pace of battle is defined to be: ¥@b. When our discrete model is fit to 


an entire battle trajectory the pace of battle provides an overall measure of the 
violence with which the battle was fought. It may be of even greater value when 
calculated at each individual time increment [n,n+ 1] of the battle thereby providing 


a picture of how the intensity of battle varied over time. 
2. Relative Combat Power 


The relative fire power is defined to be: When our model is fit to an 


& | o, 


entire battle trajectory, calculation of the relative combat power gives the average 
relative combat power of the two opposing forces. Calculation of the relative combat 
power during each time increment [n, n+ 1] again provides a picture of how relative 
combat power varies over time. This can be viewed as how the tide of battle 


fluctuates with time. 
3. Probability of Victory 


Based on our model (equations 28), the probability of victory for the X force 


y2 
can be approximated by: P[X Wins] =~ P = < . Note that this is only an 
X, 


0 


approximation since we are ignoring the noise terms in equation 28. This measure 
allows us to measure how the distribution of the relative combat power and initial 


force ratios impact the battle outcome. 
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4. Elasticity 
Our final measure, although not coming directly from our model, is composed 


of the same data elements as the previous measures. It is known as the elasticity 


(Ref. 9]: 











Ke 
ee 

"AY, AY, 
Y 


Calculating elasticity at each time interval, we obtain a picture of which force is 
winning or which force has the momentum. Elasticity is similar to relative combat 
power in that it also provides a picture of how the tide of battle varies over time. 
5. Best Measures 

Of these four measures, using a temporal trace of the pace of battle in 
conjunction with the relative combat power provides the most descriptive graphical 
representation of each battle process; providing a view of how the violence and tide 
of battle plays out over time. The probability of victory appears to provide us with 
a useful measure for analytical comparison. 

Considering possible ways of computing the pace of battle and the relative 
combat power suggests two methods of comparison. 

1. If we calculate each measure at each time interval then we get a 

temporal trace of the pace and relative combat power during the battle (see Figure 
2). To compare these traces for two different battle trajectories we measure the 


distance between them (see Barr, Weir, Hoffman [Ref. 10]). 
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Figure 2 Temporal trace of pace and relative combat power. 
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Le If we calculate the average relative combat power 


N-1 


B, 
a : la i 
A N 
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where i = index number of the battle trajectory, 
N = number of battle trajectories in the battle data set, 
over each battle we might be able to compare these measures directly using the 
Standard t-test since the distribution of relative combat power appears empirically to 
be normally distributed. 
A methodology for using the probability of victory is described in detail in the 


next section. 


G. COMPARING THE HIGH RESOLUTION MODELS 
1. Methods of Comparison 


Regardless of the fitting technique used for estimating the parameters of the 
model (see Chapter III), the process of fitting the proposed analytical models to the 
simulation data results in a set of estimates of the realizations of the coefficient 
random vectors. We desire to use the discrete model with the estimated coefficients 
to compare the high resolution simulation models. To do this, we need to identify 
measures that can be used to detect differences or similarities in the underlying 
structure of the high resolution models. Direct comparison of the coefficient random 
vectors can rely on independence between the coefficients and on sufficiently large 
data sets. To avoid having to make an undesirable independence assumption, and 
to utilize smaller, less costly data sets, we examined two indirect approaches for 
comparing our simulation models. 

(1) Fighting simulated battles. This involves choosing coefficients from 


the empirical distributions of the fitted coefficients, and using them to run numerous 


Ze 


replications of a simulated battle. The MOEs from these simulated battles are then 
compared. 

(2) Comparing Probability Surfaces. Using the estimated coefficient’s 
empirical distribution we generate probability-of-win surfaces for various initial force 
ratios. Surfaces generated from two different battle data sets are compared to each 
other numerically to obtain a difference measure. A randomization test is then used 


to compare the two battle data sets. 
2. Fighting Simulated Battles 


Our first method uses the fitted analytical models to generate larger data 
samples for comparing the MOEs from each model being considered. 

We draw pseudo coefficients from empirical distributions of the estimated 
coefficient random variable realizations. These coefficients are then substituted into 
the surrogate analytical models and a single battle is fought using this model, yielding 
a battle trajectory. We then draw a second set of coefficients and repeat the process. 
This continues for a predetermined number of replications. 

From each simulated battle trajectory generated in this fashion we observe 
the MOE of interest. The empirical distributions of the MOEs corresponding to the 
models of the two simulation models are compared directly. 

One difficulty with this method is determining the number of replications 
required to distinguish between two different battle data sets from two different 
simulation models while still being able to recognize when two battle data sets are 


from the same simulation model. We attack this problem by computing the sample 


means (fi, and fi,), and sample variances (6,7 and 6,*) of the MOE estimated 


directly with the data from the simulation model. This information is then used to 


estimate the number of observations required to obtain a confidence interval of 


jy) 


length L on the sample means as follows: 


ta 


2 , — 
2| where Gy,,, = max{6,,6,}, 


. the maximum sample standard 
ni = a Cay 


deviation of the MOE calculated with data from the two models (assuming normality 


for the distribution of the MOE), and ‘1-4 is the 1 a quantile of a t distribution 


with degrees of freedom associated with 67 


This method of comparison tended to transform (or shift) the mean value of 
the MOE upward, and to introduce additional noise. That is, the variance of the 
MOE from the simulated (fitted) system was greater than the variance of the MOE 
in the original sample. Since our intent was to tighten the confidence interval on our 
estimated MOE, the introduction of additional noise was counterproductive, making 
the comparison less sensitive. Thus, although somewhat effective for identifying 
similarities and differences in the underlying battle structure, it was determined that 


this approach should not be pursued further. 
3. Comparing Probability Surfaces 


From our analysis in Section D, we know that if (Y))’/(X,)* < b/a then the 
X force will be victorious. Since we are assuming that A and B are random 
variables, then R = B/A will also be a random variable. We can use observations 


of R obtained from our fitting process to build an empirical cumulative distribution 


function (cdf) # for R. This empirical cdf determines a probability surface Fp over 


the domain {(Xp,Yo) 2 X;SXoSX, » WiSVoSYy } 
2 
peu 
Yo 


> 


where: F,(x,,y,) = P 








Xo = the initial fire power of the X force, 


X, = the lower bound on the initial X fire power values, 
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x, = the upper bound on the initial X fire power values, 
and similarly for the Y values. 
Suppose we wish to compare two sets V and W of battle data (e.g. one 
consisting of several battle trajectories from Janus(A) and one consisting of several 
- NTC battle trajectories). A natural measure for this comparison is the "volume" 


between the two surfaces: 


x 


Yu u 
v(V,W) = > > | Fy (Xo) - Fy Go) | Ax Ay, 
Yo= Yj % = % 
where Ax, Ay define the resolution of the partition of the domain of Fp (Ax = Ay = 


1 was used in our trials). 

A question is, "What does a given value of v(V,W) mean?" For example, 
does a value of 800 indicate that the battle sets are the same or different? To 
address this question we first compare probability surfaces with a randomization test 
aS follows. 

Given our surrogate model (equation 28) and two battle data sets V and W 
consisting of n battle trajectories each: V = {v,, v5, ..., v,} and W = {w,, Wo, ..., Wat, 
where v, = the i'" battle trajectory from simulation model V, and w, = the i" battle 


trajectory from simulation model W, we estimate the coefficients a and b for each 


battle trajectory in V and W. We denote these estimates bya,, a oe and Oe 


3 
Ww; 


where 4 _ = the estimate of a obtained by fitting battle trajectory v,, and the other 


v 


estimates are defined in the same manner. 


We define the sets R, and R, as follows: 


a 

4 

R, = " r : Lo ' where LS ra 

Vi 

and 

G,, 
R= ir ,f,*'°, 7 \; where r, = — 
Ww w,? Ww,” > w, > Ww; rs - 


Vi 


2 


We now define an empirical cumulative distribution function F, for a set of 
observed values of the random variable R; {Fay Fay os rn} where 
ray Slay S °° Slay 


max {i | ras r} 


FQ) = [Ref. 13]. 


Examining the xp-y) plane determined by a finite range of values for the initial fire 


power of the X force (X,) and Y force (Y,), that is, the plane 
{ (<ps¥o) |xs2, 5x, ¥, Soe } where x,, X,, Y), Yu are integer lower and upper bounds 
on the range of values for x, and yo, we partitioned this plane into a grid of square 
intervals 

{ @y) |x, +k <x < x, +(Kel);, y, ey y,+(k+1)} ,k = 0, 1,... , m-1; wheremige 
Xi - X= Yu Yi: 


We now use the definition of F, and the partition on the domain of F, to 


define a probability-of-win surface $. Given a set of observations of the random 


A 


variable R, {r,, ras ctea, ry} of the ratios } choose a random subset R' of R. 


G 


We then define 


2 


Z 
a ie 
Sp: = Or Vor ee SX SX; YI Sm he Fel 3 
‘(0 


We define our test statistic to be: 





As xy ye yo 
rN 0 A 0 

v(x Su) = ES Vg (| Ae [22] 
Yory; X= *; Xo Xo 
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V (Sars Sz) is the volume between the probability-of-win surfaces generated by the 


sets R, and R, » where R; and R; are sets of observed values of the random variable 


R. 


To determine if the two simulation models V and W are different, we first 
compute Tf = V(SR ) 5 Ry) This is the observed value of our test statistic. We then 
form R = R, U R,,. Let P be the set of all possible partitions R into two subsets R, 
and R, such that both R, and R, contain n elements and R,  R, = o. A specific 
partition p € P can be built by selecting n elements of R at random (without 
replacement) as the set R, and letting the set R, be R\R,. We can then compute 


V (Sp,» 5p). Our randomization distribution is defined to be the set of values 


v (Sr,> N R,) for all p e P. We compare our observed value T with a sample from the 


randomization distribution, obtained by randomly choosing M partitions p e€ P and 


evaluating v (S RY S R,) for each p. If the value of T is in the right tail of this sampled 


distribution (i.e. an extreme value) we infer that V and W are from different 
simulation models. And if T is in the left tail of the randomization distribution we 
infer that V and W are from the same simulation model. 


Since the number of values in the randomization distribution may be as high 


as i a = Gn it would be impractical to compute each value. However, it is 


7) 4(n!) 
reasonable to take a sample from this distribution for testing the observed value of 
T [Ref. 1]. 
In our studies, this technique proved to be effective (results are given in 
Chapter IV) at identifying difference in battle data sets V and W, when V and W 
came from different generating systems (see Chapter IV , Section A for a description 


of the generating systems used), while still recognizing similarity, when V and W 


ae | 


were generated by the same system. This technique is demonstrated in Chapter IV. 
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II. FITTING TECHNIQUES 


A. INTRODUCTION 

Our approach to finding an appropriate technique for fitting the coefficients of 
a system of simultaneous difference equations has been experimental. We first 
examined the simulated annealing technique suggested by Ingber [Ref. 2]. We then 
considered methods that required less computation time. This chapter is a 


chronological account of our research, we thus present our results in that order. 


B. SIMULATED ANNEALING 

The simulated annealing algorithm was developed by Metropolis in 1953 [Ref. 
11] to simulate the physical annealing process studied in statistical mechanics. It was 
first suggested as a technique for solving combinatoric optimization problems in 1983 
by Kirkpatrick, Gelatt, and Vecchi [Ref. 12]. The simulated annealing algorithm 
combines local optimization techniques with a random walk process to reduce the 
chance of becoming trapped in a local optimum. The algorithm shown in Figure 3 
begins with an initial solution and generates a proposed neighboring solution at 
random. If the proposed solution is better than the current solution (i.e. a downhill 
move) then the proposed solution is accepted as the new current solution with 
probability one. If the proposed solution is worse than the current solution (i.e. a 
uphill move) then it 1s accepted with a probability based on the magnitude of the 
uphill move and the current value of the control parameter (referred to as the 
temperature of the process). Thus, small uphill moves are more likely to be accepted 
than large ones, and all uphill moves are more likely to be accepted at higher 
temperatures than at lower temperatures. The process is run for a large number of 
repetitions at each temperature, and reductions in temperature are made according 


to some "cooling schedule." Changes to the initial solution, initial temperature, and 
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cooling schedule can have dramatic effects on the convergence properties of the 
algorithm and most often must be arrived at by experimentation. 

Simulated annealing’s main utility is in solving hard combinatorial optimization 
problems for which: 


e no good problem specific heuristic algorithms have been developed (simulated 
annealing does not compete well with problem specific techniques, like those 
developed for the traveling salesman problem, see Johnson, Aragon, McGeoch, 
and Schevon (1989) (Ref. 14]); 


. Get an initial solution S. 
. Get an initial temperature T>0. 
. While not yet frozen do the following. 
3.1 Perform the following loop L time. 
3.1.1 Pick a random neighbor S?’ of S. 
3.1.2 Let A = cost(S’) - cost(S). 
3.1.3 If A < 0 (downhill move), 


Seley = Se 
3.1.4 If A > 0 (uphill move), 
Set S = S’ with probability e4”. 
3.2 Set T = rT (reduce temperature). 
4. Return S. 


Source: JOHNSON ET AL. [Ref 13]. 





Figure 3. Generic Simulated Annealing Algorithm. 


e there is limited application for the formulated problem (if the problem 
formulation has wide applicability efforts toward developing an effective 
problem specific algorithm should be more beneficial); 


e the solution space is well understood (since the annealing algorithm is highly 
parameter dependent, a strong understanding of the solution space is critical 
for ensuring convergence of the algorithm to a global optimum). 


In a proposal to the U.S. Army TRADOC Analysis Command [Ref. 3], Ingber 
proposed that a variation of the conventional simulated annealing algorithm, referred 


to as Very Fast Simulated Reannealing, could be used effectively to fit systems of 
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differential equations to battle data from various high resolution simulations and 
NTC. Our initial experiments with simulated annealing were designed to use a 
somewhat simplified version of Ingber’s algorithm to assess his claim. We first 
attempted to fit individual battle trajectories and various subsets of battle data sets 
to a system of differential equation of the following form (our difference equation 


model was developed later): 





au = AptankpranKBEANK + OpranxzapcBAPC + GeranxsrowBTOW 
anes = Gespcarane BANK + Gpspc pape BAPC + Gpanc prow BTOW 
aos = OpranxpTankREANK + Garang papcRAPC 

E — © = Gpapcprane RTANK + Opyoc pap RAPC 

am = Gprowrrane RTANK + Gprow pape RAPC 


where RTANK denotes the number of red tanks involved in the battle at time t, 
ARTANKBTANK IS the attrition rate of red tanks by blue tanks, and similarly for the 
other variables and coefficients (note here that APC denotes an armored personnel 
carrier, and TOW denotes the Tube-Launched, Optical-tracked, Wire-guided antitank 
weapon system of the U.S. Army). 

After numerous experimental runs of the simulated annealing algorithm it 
became very clear that regardless of the parameters used in the algorithm, when a 
small number of battles were used, the results were quite unstable. In particular, 
Starting with one seed for the pseudo-random number generator we would get one 


solution and starting with a different seed we would arrive at a different solution 
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using the same data. When a larger number of battles were used the computer time 
required became unreasonable (on the order of 24 hours on a 386-based PC running 
at 33 Mhz when fitting ten battle trajectories). | Based on these results we 
temporarily abandoned simulated annealing and searched for quicker and more 
Stable techniques. We discuss several such techniques in the remainder of this 
chapter. We found during our experiments with various fitting methods that it was 
reasonable to specify relative fire power indices for each weapon system, thus 


allowing us to fit a simpler system of differential equations given below: 


oe =O Lae 

dt 

LOO. 

dt e 
where, 


X = feraygRTANK + fespcRAPC, 


Y = forygBTANK + fospcBAPC + feroyBTOW. 


Experimental runs using simulated annealing for fitting this system of equations 
were both stable and much quicker than those experienced previously. However, we 
observed that the simulated annealing process simply converged to the regression fits 
discussed in Section E, and took about 10000 times longer to do so. Thus, simulated 
annealing was found to be an unsuitable technique for fitting systems of equations 


of the form in which we were interested. 


C. STEEPEST DESCENT 

While trying to determine why simulated annealing was not working well, we 
made two observations regarding our model. First, Lanchester models assume that 
every weapon system on the battlefield can see and engage every other system on the 


battlefield. This was not true of either the battles observed at NTC or those 


a2 


generated by Janus(A). Second, since the models under consideration were attrition 
models, their solution spaces were convex, and thus a least squares measure (i.e. the 
sum of the square differences) of the fit of the data to these models would also 
induce a convex space. Thus a simpler fitting technique such as a standard steepest 
descent algorithm should be effective. 

At this point we consulted with analysts at the TRADOC Analysis Command 
(TRAC) and determined that although battle data showing which systems could see 
and engage which other systems was not currently available, the Janus(A) system 
could be modified to produce such data. Thus, we proceeded by generating 
simulated battle trajectories from mathematical models to use in testing our fitting 
techniques (see Chapter IV, Section A). 

We began testing the straight-forward steepest descent algorithm shown in Figure 
4, using this simulated data (FORTRAN code is in Appendix C). 

This algorithm proved to be stable with respect to the quality of solution, but not 
with respect to time of solution. That is, regardless of the starting solution , the 
algorithm would converge to the same solution for a given accuracy level; however, 
the time required to do so varied from one minute up to four hours. Although this 
algorithm was far more useful than the simulated annealing approach, the idea of 
using non-linear fitting techniques led us to try using commercial non-linear 


programming software to fit our system of equations. 


D. NON-LINEAR OPTIMIZATION 

The convex nature of our problem encouraged us to believe that we might be 
able to use commercial non-linear solvers to fit our equations to the battle data. We 
thus formulated our problem on the General Algebraic Modeling System (GAMS) 
using the MINOS 5.1 solver (see code in Appendix D). We made experimental trials 


using both a least squares measure (i.e Euclidean distance) and a Chebychev (or 


aR) 


minimax) criterion. 


Least Squares Objective Function: Min }~ (6, + 6 a 


Chebychev Objective Function: Min ( max{ 5.1 + 15, 1}) 
Subject to: 6, = (%,,, - x) - @y, + Z) 


xX, i 


oF > Opgui= yeast zy) 


a 


Solve for: Gop bopr Cx)opr? ©) opr 


. Get an initial solution S. 
. Determine an initial grid definition size. 
. While not yet close enough: 
3.1 Compute cost for each neighboring point on the grid. 


3.2 Set S’= neighboring solution with min cost. 
3.3. IF IIS - S’l < accuracy THEN refine the grid size. 
. Return S = S’, 





Figure 4. Steepest Descent Algorithm. 


where 6, = difference between the predicted attrition and the actual attrition. 


Our results were encouraging. The solutions were stable and were found 
consistently in just under one minute for single battles. We also observed graphically 


that the least squares criterion provided a smaller variance on the fitted coefficients. 


E. LINEAR REGRESSION 

Our GAMS formulation with a least squares objective function was essentially 
providing a simultaneous linear regression fit of the battle trajectories to the 
equations in the analytical model being fit. This led us to ask when this would be 


the equivalent or nearly equivalent to doing linear regression on each of the 


34 


individual equations separately. The answer is that if we assume the coefficients in 
the equations are independent, then we could simply fit the equations separately 
using regression. 

Considering the method (Chapter IV, Section A) used to generate our simulated 
battle trajectories (which held the attrition and noise coefficients constant for the 
duration of each individual battle), this independence relationship was true for the 
individual battle trajectories (this would not be the case for the data in general). We 
did regression fits on the individual equations and found the two methods equivalent 
(i.e. the fitted parameters were nearly equal) up to a factor that could be explained 
by noise induced by the generation of integer data points. 

Thus, if we are willing to make the assumption that the equations are 
independent for individual battles, we reduce our fitting time from about one minute 
per battle to about 0.1 seconds per battle. Whether this makes sense to do on actual 


battle trajectories needs to be examined once the data become available. 


F, CONCLUSION 

Our experiments with the non-linear optimization approach of using a GAMS 
model with a least squares objective function indicate that this approach is the most 
effective way of fitting a system of difference equations to battle trajectories, 
providing consistent solutions in a timely manner. However, experiments should be 
made to consider whether assuming independence between the coefficients in each 
equation of the surrogate model is justified, thus allowing us to use the faster 
independent regression technique discussed in the last section. We use the 


independent regression technique for an example in the next chapter. 
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IV. EXAMPLE OF METHOD 


A. THE BATTLE DATA TRAJECTORIES 

To test the methodology, we first generated battle trajectory data based on a 
discrete version of Lanchester’s aimed fire model with additive noise term (equation 
29 below). Three different systems of equations (all of this same form) were used 
for this purpose. These generating systems differed only in the fire power indices 
assigned to each weapon system. Each system was used to generate two separate 


collections of 20 battles each. The form of the generating equations is: 


Oe a X, a = Axe i Z,. 
Y,4-Y, = -BX,+Z,, (29) 


where 
X, = Traine RI ANKS + eure RAPC,, 


Y. = fore BIANK, + fespc BAPC, + feroy BTOW, - 


and A and B are constant throughout a single battle trajectory. Here, fpranx is the 
fire power index for the red tanks, and similarly for the fire power indices of the 
remaining weapon systems (values shown in Table 2 below). RTANK, is the number 
of red tank systems active in the battle at the end of time increment n. RTANK, is 
computed as follows from X,. The other weapon system levels (RAPC, BTANK, 
BAPC, and BTOW) are all computed in the a similar manner. At the end of each 
time interval [n, n+1], the attrition sustained by the X force is divided among the 


two X weapon systems as follows: 


RTANK,,, = RTANK, + (X 


n+l 3) 


tfpwry - Srranx 
(tfpwr,)(Xsys - 1)’ 





RAPC,,, = RAPC, + (X 


n+l ) 


pwry - Srarc 
(fpr) (Xsys - 1) 
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where 
tfpwry = feranx * Srarc 
Xsys = the number of different weapon system in X force, 
and similarly for the Y force. We then round RTANK,,, and RAPC,,, down to the 
next integer value. 
Also, [A Z, B Z,] is distributed multivariate normal with variance-covariance 
matrix: 
0.00004 -0.00200 -0.00001 0.00010 
-0.00200 0.50000 -0.00060 0.40000 
~ |-0.00001 -0.00060 0.00008 -0.00600]’ 
0.00010 0.40000 -0.00600 1.00000 | 


and mean pz = [ 0.066 0.000 0.055 0.000 J. 


Initial force levels were 


RTANK(0) = 80 BTANK(0) = 60 
RAPC(0) = 160 BAPC(0) = 120 
BTOW(0) = 40, 


The fire power indices f were set as shown in Table 2: 


TABLE 2. FIRE POWER INDICES. 
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A realization of [A, Z,, B, Z,] was chosen and a battle trajectory consisting of 12, 
five minute intervals was run out using the generating equations (29). Another 
realization was then chosen and another battle trajectory was generated. In this 


manner 20 battle trajectories were generated for each battle data set. A sample 


battle trajectory generated with these equations is shown in Table 3. 
TABLE 3. BATTLE TRAJECTORY SAMPLE. 


Battle Trajectory (Battle 1 Set A) 
BTANK BTOW 
60 40 





B. SELECT MATHEMATICAL SURROGATE SYSTEM 

We had decided to use the Lanchester aimed fire model in our initial tests. 
However, we had to decide whether to include the additive (random) noise term. 
The concern was to choose the model that would best represent the underlying 
structure of the data represented by the generating system. The best candidate was 
determined by running a small sample test. Five battles were selected at random 
from battle set A (see Table 3) and were fit to the discrete form of the Lanchester 
aimed fire model with and without the constant noise terms. The optimization 


model (procedure) used for estimating the parameters in the following surrogate 
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model without additive noise terms 
Xie = Xx, = |-4 4 e,{n)|Y, » 


¥,..-¥, =(-b + e(n)]X,, (30) 


is given below. 
Proportional Procedure (no additive noise term): 


Minimize: }° (6, 


2 2 

+ 6, 
subject to: 6, 2 (X,,, - X) - @Y,) 
8, 2 (¥,.. - ¥) - GX) 


Solve for: Gopp b a 
where 5y and 6, are as in Chapter II. 


Additionally, the surrogate model with the additive noise term (equation 27) was 
fit using a one-stage and two-stage fitting procedure. The one-stage procedure was 
the non-linear programming technique describe in Chapter III. 

One-stage Procedure: 
once ° i 2 2 
Minimize: d (6, + 8) 


subject to: 6, 2 (X 


int ~— A) - (AY, + 2) 
8, 2 (Vi - ¥) - OX, +2,) 


A 


Solve for: a opr Dopp qx pT’ <y opr 


The two-stage procedure used the same basic technique; however, in stage one the 


noise terms were fixed at zero and the attrition coefficients were fit (this stage is the 
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same as the proportional procedure). In the second stage the attrition coefficients 


were fixed at the optimal value found in stage one and the noise term was fit. 


Two-stage procedure: 
Stage One Problem: 


Minimize: pre, + 6, | 
subject to: 6, 2 (X,,, - X) - @yY) 
6, 21, — (oe) 


Solve for: Gopp Oopr 


Stage Two Problem: 
Minimize: }° (6, 


subject to: 6, 2 (X,,, - X) -((@5,,¥ + z) 


oi 


~~ 


oF = (Y., - ¥) - GoprX; + 2) 


Solve for: Z a 
uk <x opr’ <y opr 


The results are shown in Tables 4, 5, and 6, where cost is taken to be the square root 


of sum of square differences of the predicted and observed attrition from the model 
(equation 29). 
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TABLE 4. ONE-STAGE FITTING. 


Results of One-Stage Fitting Procedure 


re EC 





TABLE 5. TWO-STAGE FITTING. 


Results of Two-Stage Fitting Procedure 
BATTLE 





TABLE 6. PROPORTIONAL FITTING. 


Results of Proportional Fitting Procedure 
BATTLE 
5378 ig] saa? | 10991 


0.061 0.049 0.048 0.054 0.056 
0.042 0.048 0.037 0.049 0.049 


Examining the costs (sum of squares error) in the tables above we concluded that 







the one-stage fitting procedures was the best (minimum error) of the three examined. 
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We were however interested in obtaining a fitted model that best represents the 
underlying structure of the data that we were fitting. Thus, we also looked at the 
Euclidean-norm of the difference between the fitted coefficients and those used to 
generate the data. Let a = [A, Z,, B, Z,); define & to be the parameters estimated 
by the fitting procedure (take Z, = Z, = 0 in the proportional procedure). Let 


a, = c ar Hz» Bes Hz |, the mean values of the parameters used in the generating 


equations to simulate the battle data _ sets. We then compute 


+ (Z, 





la - @ = i (6 2 Ls) + (2, a bz): These values are 


(4 - Hy) bz) 


Shown in Table 7. 


TABLE 7. NORM OF DIFFERENCES IN COEFFICIENT VECTORS. 


Euclidean-norm of difference in Coefficients 
A2 EE: 


refs 1.78 


We also examined the battle trajectories for each of these five fitted surrogate 





models graphically, by plotting the force levels versus time for the fitted system and 
the generating system (see Figures 5-7 for an example). 

As a result of the analyses we concluded that although the one-step fit was 
providing the best fit in the least squares sense; when we examined the results using 
the proportional fit on each individual equation (in equation 28) independently 
(Figure 8) the fit to the underlying structure (i.e. the generating equations) was 
providing a graphical representation almost as good as the one-stage procedure. 
Since the proportional procedure used was much quicker than the one-stage 


procedure we decided to use it. The following optimization model is the procedure 
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Figure 5. Battle Trajectory Plots. (1-Stage vs Gen) 
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Figure 6. Battle Trajectory Plots. (2-Stage vs. Gen) 


44 


Time Sequence Pilot ——  TGEN RED 
——  TRPRO.RED 


(Poy @) 


20 





90 
0 2 4 6 8 10 Te 
Time 
Time Sequence Plot SS RGeNy Bele 


--— TRPRO.BLUE 





eo 
210 
isle 


a7 0) 


gio 


Figure 7. Battle Trajectory Plots. (Pro vs. Gen) 
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Figure 8. Battle Trajectory Plots. (Pro-I vs. Gen) 
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used. 


Proportional independent procedure (PRO-I): 


NoWex aso. 
i+] i 
py 


i=0 


Aes 
N 
ee a = 
A AE X. 
en Nat es) 
N 


We therefore used the discrete difference equation model without additive noise 


(equation 30) for the remainder of the analysis. 


C. FITTING THE MODEL TO THE DATA 
Based on the analysis in the preceding section, the model was fit to the battle 
data sets obtained from the generating systems using the proportional independent 


method. The coefficients fit to the first five replications of each data set are shown 
in Table 8 below. 
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TABLE 8. FITTED COEFFICIENTS. 


Sample of Fitted Coefficients 


iE 0.0586 0.0502 0.0493 0.0406 0.0499 0.0414 | 
0.0426 0.0505 0.0499 0.0582 0.0417 0.0466 

Za 0.0469 0.0489 0.0383 0.0379 0.0387 0.0394 
0.0487 0.0530 0.0548 0.0592 0.0451 0.0474 


So 0.0520 0.0533 0.0409 0.0447 0.0428 0.0449 


0.0400 0.0396 0.0515 0.0459 0.0411 0.0391 

0.0509 0.0573 0.0422 0.0469 0.0423 0.0479 

0.0524 0.0499 0.0580 0.0572 0.0469 0.0464 

[4 0.0600 0.0572 0.0505 0.0463 0.0512 0.0476 
: 0.0503 0.0470 0.0584 0.0551 0.0467 0.0446 
0.0471 0.0462 0.0479 

0.0555 0.0448 0.0452 


D. EMPIRICAL DISTRIBUTION FUNCTIONS 
Graphical analysis of the fitted coefficients from the battle data sets (of 20 battle 





trajectories each) show no clear fit to a known probability distribution, thus we 


turned to the empirical distributions. The empirical cumulative distribution functions 


for the ratio of the estimated coefficients Lu where determined for each battle (see 
a 


Figure 9 below). 

Using the procedure outlined in Chapter II we generated the probability surfaces 
numerically over square intervals and calculated the difference between surface 
heights at the corner of each square interval {X < x) < X+1, Y < yy) < Y+1}. The 


observed values of our test statistic T (volume between surfaces compared) of the 
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EMPIRICAL CUMULATIVE 01STRIBUTION FUNCT ION 


Probab! I Ity 
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Relative Combat Power 


Figure 9. Empirical CDF for B/A (Battle Data Set A). 
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battle set comparisons are listed in Table 9. The notched box plots in Figure 10 
provide a graphical portrayal of samples of size 50 of the randomization distribution 
of these differences and the location of the observed test statistic within this 
distribution. Examination of these box plots clearly delineates the differences 
between the dissimilar battle data sets, and the similarity of the equivalent battle 


data sets (i.e. A and B, C and D, and E and F as defined in Table 3). 


TABLE 9. OBSERVED VALUE OF THE TEST STATISTIC T. 


——_ 2s LLL 
| 2940 | 2849 | 1079 | 98 


1769 Webay 





E. ANALYTICAL COMPARISON OF BATTLE DATA SETS 

Having calculated the differences between probability surfaces generated from 
the empirical distribution functions of random samples of the ratios of fitted 
coefficients, we now use the 50 samples obtained from the randomization distribution 
of the surfaces differences and the observed value of our test statistic (T) to compare 
the battle data sets from the two simulation models being compared. To do so we 
simply observe the relative location of the observed statistic in the randomization 
distribution by calculating the percentage of realization values greater then the 
observed value. This is equivalent to determining the percentile in which it lies. If 
the percentage is large then we do not reject a hypothesis that the two battle data 
sets came from similar (possibly the same) simulation models. However, if the 
percentage is small then we reject that hypothesis. Since our trials resulted in 


percentages well out in the tails of the randomization distributions (see Table 10), 
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we did not postulate any conclusions with respect to an exact cut off percentage for 


rejection. 


TABLE 10. ANALYTICAL COMPARISONS OF RESULTS (THE NUMBERS 
INDICATE THE PERCENTAGE OF SAMPLE VALUES OF THE 
RANDOMIZATION DISTRIBUTION GREATER THAN THE TEST 
STATISTIC). 


Analytical Comparison 


DO NOT | REJECT | REJECT REJECT | REJECT 

REJECT 

0.96 0.00 0.00 0.26 0.14 
REJECT | REJECT REJECT | REJECT 
0.00 0.00 0.08 0.08 


DO NOT | REJECT | REJECT 
W(C, ) REJECT 
0.74 0.10 0.04 
REJECT | REJECT 
v(D,_) 
0.04 0.02 
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SAMPLE REFERENCE O!STRIBUTIONS OF VOLUME BETWEEN SURFACES 


(Numbers indicate observed value of the test statistic) 
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Battle Set _ vs Battle set _ 


Figure 10. Box Plot Comparison of Surface Differences. 
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V. CONCLUSIONS AND RECOMMENDATIONS 
A. INTRODUCTION 


This thesis has presented a method for comparing highly complex combat models 
using a simple analytical surrogate model. In doing so, we have evaluated a number 
of fitting techniques and other comparison techniques, and we have developed a 


discrete analog to Lanchester’s aimed fire model. 


B. SURROGATE MODELS 

We have concentrated here on Lanchester-like models of combat; however, other 
models could be studied. Non-homogeneous models which include area fire weapons 
(such as artillery and mortars) were not studied in this thesis but should be studied 
if this methodology is to be useful as an analytical tool. Models specifically 


addressing human factors should also be explored. 


C. FITTING TECHNIQUES 

Although we found that the simultaneous least squares (GAMS model) technique 
most effectively represented the underlying structure of battle data sets fitted to our 
model, all of the fitting techniques considered have utility for certain types of 
problems. Simulated annealing is a very time consuming procedure, but could be 
useful if the simulation models being studied allow reinforcement. A _ useful 
algorithm for fitting non-linear, non-convex continuous functions using simulated 
annealing is discussed in Brooks and Verdini [Ref. 15]. 

GAMS provided a very flexible experimental tool for fitting our models. 
Numerous possible fitting techniques can be written in GAMS in an efficient manner. 

We expect that attempts to use multiple regression (as in our proportional 
independent procedure) with more complex models will experience difficulty with the 
independence assumption. It has, however, provided a powerful tool for fitting a 


large number of battles in a very short time. 


Bis, 


D. RECOMMENDATION 

Comparing high resolution models of combat continues to be a challenging 
problem. We have seen here that it is feasible to use very simple analytic surrogate 
models to compare machine generated battle data sets. There is also a lot of work 
to be done in applying this methodology to actual data sets generated by Janus(A) 
and at NTC, and in investigating various surrogate models. Successful comparison 
of these two highly complex models of combat should be of great benefit to the Army 


analysis community. 
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APPENDIX A 
MATRIX SOLUTION TO DISCRETE MODEL 
In Chapter III we developed the 2 x 2 case of the discrete difference equation 
model. In doing so we combined the fire power from the various weapon systems 
into one term by using the fire power indices. Suppose instead we wish to consider 
a model where the attrition of each weapon system and its contribution to the 
attrition of all other weapon systems are considered separately. The discrete model 


for this system would be given by: 


Nynta > Ank = -4y Yin -Ai Yor ® © © -Aim Yme > 
Xo x41 7 Mon = “Ag, Yan “B22 You © © © -Bam Yank » 
@ 
@ 
@ 
eet te Sat Lik Snr Yon © 8° aam Yk » 
Lk+1 7 Yur = -b,, Xi k -by» Xo mates Din k » 
eel Yo = -b,, Xi “by, Xk eee -d,, nk » 


Mee e Rann = Oi Xi, Dn? X25 eee -b Xa : 


mn 
where, 
X;, = the number of X force weapon type 1, active in the 
battle at the end of time increment k, 
Y,, = the number of Y force weapon type Jj, active in the 
battle at the end of time increment k, 
a,; = the rate at which Y force weapon type J kills 
X force weapon type 1, 
b,; = the rate at which X force weapon type i kills 
Y force weapon type j. 


>) 


Or in matrix notation: 


10C0QG0::-:. Qg Ub, es = 





1,k+1 lm 
Xen O10 ++: 0 -ai2 +--+ -a, 
eel 0 OU =a eee 
ea | | ee -b, 100---0 
Yee -b21 +--+ -b 0 1 0---0 
Dekel | bt ae 0 0 1 





In block notation we can write: 





Yoel ; Bot Y, | 

















OT 


Feat 


IA 
BI 





= RF,_; where R = 


The solution to this system is: 


F(k} = R* F(O). 


Associated with R are its eigenvalues {,, ,, ¢ © © , A,4m} and their corresponding 


eigenvectors {A,, A,,e e ©, A,,_}. These eigenvectors form a basis for R"*™ so we 
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1k 


2,k 


m,k 


can then write F(Q) as a linear combination of these eigenvectors: 
FO) =c,A, + 6,4, +--+ +e 


n+m~"n+m° 


So finally, the solution to the discrete system can be written: 


F(k) = (5) (A,)* A, we (A,)* A, eee Cae 
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APPENDIX B 
SIMULATED ANNEALING CODE 
PROGRAM siman 


C** Simulated Annealing Algorithm implemented for the 
C** discrete analog to Lanchester’s Aimed Fire Model. 
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REAL  cred(100),cblue(100),red(100),blue(100),c(4), 


& acost,ccost,fnacc,fniter,pacc,temp,cool,minpct 


INTEGER n,maxcyc 
COMMON /vars/cred,cblue,red,blue,n 


INTEGER nacc,niter 
REAL tcost,mesh(4),a(4),at(4) 
DOUBLE PRECISION dseed 


dseed = 18564245.0d0 


DO31i=1,4 
a(i) = 0.0 
c(i) = 0.0 
CONTINUE 


mesh(1) = 0.0001 
mesh(2) = 0.001 
mesh(3) = 0.0001 
mesh(4) = 0.001 


temp = 100.0 
cool = 0.95 

maxcyc = 100 
minpct = 0.01 


CALL readit 


CALL evalfn(a,acost) 
ccost = acost 
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C** Begin Simulated Annealing Algorithm 
DO 36k = 1, maxcyc 
nacc = 0 
niter = 5000 
DO 35 j = 1, niter 
CALL nbhd(a,at,mesh,dseed) 
CALL evalfn(at,tcost) 
IF (tcost .le. acost) THEN 
DO 331=1,4 
a(i) = at(i) 
IF (tcost .It. ccost) c(i) = a(i) 
33 CONTINUE 
acost = tcost 
IF (acost .It. ccost) ccost = acost 
nacc = nacc + 1 
ELSE 
CALL unif(dseed,u) 
IF (u .gt. exp(-temp*(tcost-acost))) THEN 
DO 341=1,4 
a(i) = at(i) 
34 CONTINUE 
acost = tcost 
nacc = nacc + 1 
ENDIF 
ENDIF 
35 CONTINUE 
fnacc = FLOAT(nacc) 
fniter = FLOAT(niter) 
pace = fnacc / fniter 
C** Geometric Cooling 
temp = cool * temp 
IF (pacc .Jt. minpct) THEN 
PRINT *,’End conditions met in ’,k,’ cycles’ 
STOP 
ENDIF 
36 CONTINUE 
PRINT *,’Maximum iterations reached’ 
201 FORMAT(1x,’A?,4(F12.5),/,1x,’C?,4(F12.5)) 
202 FORMAT(1x, ACOST?’,F12.5,) CCOST?,F12.5, 
& PCT ACCEPTED:’,F8.3,/) 
END 
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SUBROUTINE readit 

REAL — cred(100),cblue(100),red(100),blue(100) 
INTEGER n 

COMMON /vars/cred,cblue,red,blue,n 


INTEGER ulist,ufpwr,udata,uout,d(5),numobs,i,time 


REAL fpwr(5) 
CHARACTER*9 bDfile 
ulist = 8 

ufpwr = 9 

udata = 10 

uout = ll 


OPEN(UNIT =ulist, FILE =’slist.dat’?»,STATUS =’old’) 
OPEN(UNIT =ufpwr,FILE =’fpwr.dat’,STATUS =’old’) 
DOAgic=-1..5 
READ (ufpwr,101) fpwr(i) 
1 CONTINUE 
n= 1 
2 CONTINUE 
C** Read list of Battle Data Set Files to be used. 
READ (ulist,102,END =4) bfile 
OPEN(UNIT =udata,FILE = bfile, STATUS =’old’) 
READ (udata,103) numobs 
DO 31 = n,n+numobs-1 
READ (udata,104) time,d(1),d(2),d(3),d(4),d(S) 
red(i) = FLOAT(d(1))*fpwr(1) + 
FLOAT(d(2))*fpwr(2) 
blue(i) = FLOAT(d(3))*fpwr(3) + 
FLOAT(d(4))*fpwr(4) + 
& FLOAT(d(5))*fpwr(5) 
IF Gi .gt . 1) THEN 
cred(i-1) = red(i) - red(i-1) 
cblue(i-1) = blue(i) - blue(i-1) 
ENDIF 
3 CONTINUE 
n = n+numobs-1 
GOTO 2 
4 CONTINUE 
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101 
102 
103 


105 


o* * 


201 


OPEN(UNIT =uout, FILE =’sa.out’» STATUS =’unknown’) 
DO 51=1,n-1 
WRITE(uout,105) i,red(i),cred(i),blue(i),cblue(:) 

CONTINUE 
n = n-l 

FORMAT(F10.7) 

FORMAT(A9) 

FORMAT(IS) 

FORMAT(6I5) 

FOR MAT(I3,5(2x,F10.7)) 
RETURN 
END 


SUBROUTINE evalfn(c,ecost) 
Sum of Square Differences (Linear Lanchester) 


REAL  cred(100),cblue(100),red(100),blue(100) 
INTEGER n 
COMMON /vars/cred,cblue,red,blue,n 


REAL ecost,r,b,c(4) 


ecost = 0.0 

DO 201i =1,n 
r = (cred(i) - (c(1)*blue(i) + c(2)))**2 
b = (cblue(i) - (c(3)*red(i) + c(4)))**2 
ecost = ecost + r + b 

CONTINUE 

ecost = SORT (ecost) 

RETURN 

END 
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Crs 


Cr. 
Ce 


SUBROUTINE nbhd(b,bt,m,dseed) 
Chooses random neighboring solution. 


REAL b(4),bt(4),m(4),u,s 
DOUBLE PRECISION dseed 


CALL unif(dseed,u) 
CALL unif(dseed,s) 

j = INT((u*4.0)+0.5) 
sign = (s-0.5)/abs(s-0.5) 
bt(Gj) = bG) + sign*m() 


RETURN 
END 


SUBROUTINE UNIF(DSEED,V) 

Uniform (0,1) Pseudo Random Number Generator (Lewis & 
Orav [Ref. 16]) 

INTEGER I 

REAL U 

DOUBLE PRECISION DENOM,DSEED 


DATA DENOM/2147483647.D0/ 

DSEED = DMOD(16807.D0*DSEED,DENOM) 
U = DSEED/DENOM 

RETURN 

END 
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APPENDIX C 
STEEPEST DESCENT CODE 


PROGRAM steepd 


C** Steepest Descent Algorithm 


of 


REAL — cred(100),cblue(100),red(100),blue(100),c(4), 
& acost,ccost,fnacc,fniter,pacc,temp,cool,minpct 
INTEGER n,maxcyc 

COMMON /vars/cred,cblue,red,blue,n 


INTEGER nacc,niter,tmax,term 
REAL tcost,mesh(4),a(4),at(4),amin(4),amax(4) 


DO 311=1,4 
a(i) = 0.0 
c(i) = 0.0 
CONTINUE 


mesh(1) = 0.0001 
mesh(2) = 0.001 
mesh(3) = 0.0001 


mesh(4) = 0.001 
amin(1) = -1.0 
amax(1) = 1.0 
amin(2) = -5.0 
amax(2) = 5.0 
amin(3) = -1.0 
amax(3) = 1.0 
amin(4) = -5.0 
amax(4) = 5.0 


precis = 0.0001 


refine = 0.60 
maxcyc = 20 
maxit = 2500 
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CALL readit 

CALL evalfn(a,acost) 
pcost = acost 

ncyc = 1 

it = 1 


C lopiotewhile looper. 
C** Finds Direction of greatest improvement in the Cost Fen. 
1000 delmax = 0.0d0 

DO 1001 term = 1, 4 

a(term) = a(term) + mesh(term) 
IF (a(term) .le. amax(term)) THEN 
CALL evalfn(a,tcost) 
delta = acost - tcost 
If (delta .gt. delmax) THEN 
delmax = delta 
tmax = term 
sign = 1.0d0 
ENDIF 
ENDIF 
a(term) = a(term) - 2.0d0*mesh(term) 
IF (a(term) .ge. amin(term)) THEN 
CALL evalfn(a,tcost) 
delta = acost - tcost 
IF (delta .gt. delmax) THEN 
delmax = delta 
tmax = term 
sign = -1.0d0 
ENDIF 
ENDIF 
a(term) = a(term) + mesh(term) 
1001 CONTINUE 
C** If an improving direction is found move in that Direction 
C** on grid [mesh] square. 

IF ((delmax .gt. 0.0) .and. (it .le. maxit)) THEN 
a(tmax) = a(tmax) + sign*mesh(tmax) 
CALL evalfn(a,acost) 
it = it + 1 
GOTO 1000 

C** If change in cost during this cycle less then set Precision 
C** We are done! 
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ELSEIF ((pcost-acost) .lt. precis) THEN 
GOTO 1003 
C** If num of cycles is < maximum cycle limit = > Refine Grid 
It. maxcyc) THEN 
ncyc = ncyc +1 
pcost = acost 
DO 1002 term = 1, 4 
mesh(term) = mesh(term) * refine 
1002 CONTINUE 
it = 1 
GOTO 1000 
ENDIF 
1003 CONTINUE 
Print *,’ncyc,acost, ’,ncyc,acost,it,maxit 
Print *,’a’,a 
PRINT *,’end’ 
WRITE(*,392) >’ COST’ ,FLOAT(acost) 


390 FORMAT(1x,A6,8x,A12) 
391 FORMAT(1x,A25,2x,F10.7) 
392 FORMAT(1x,A5,10x,F12.7) 


36 CONTINUE 
END 


SUBROUTINE readit 


REAL ~ cred(100),cblue(100),red(100),blue(100) 
INTEGER n 
COMMON /vars/cred,cblue,red,blue,n 


INTEGER ulist,ufpwr,udata,uout,d(S),numobs,i,time 
REAL fpwr(5) 
CHARACTER*9 bDfile 


ulist = 8 

ufpwr = 9 
udata = 10 
uout = 11 


OPEN(UNIT=ulist, FILE =’slist.dat’",STATUS =’old’) 
OPEN(UNIT =ufpwr,FILE =’fpwr.dat’» STATUS =’old’) 


65 


ELSEIF (ncyc 


DO 11=1,5 
READ(ufpwr,101) fpwr(i) 
1 CONTINUE 
n= 1 
2 CONTINUE 
READ (ulist,102,END =4) bfile 
OPEN(UNIT = udata,FILE = bfile, STATUS =’old’) 
READ (udata,103) numobs 
DO 31 = n, n+numobs-1 
READ(udata,104) time,d(1),d(2),d(3),d(4),d(5) 
red(i) = FLOAT(d(1))*fpwr(1) + 
FLOAT(d(2))*fpwr(2) 
blue(i) = FLOAT(d(3))*fpwr(3) + 
FLOAT(d(4))*fpwr(4) + 
& FLOAT(d(5))*fpwr(5) 
IF G .gt. 1) THEN 
cred(i-1) = red(i) - red(i-1) 
cblue(i-1) = blue(:) - blue(i-1) 
ENDIF 
3 CONTINUE 
n = n+numobs-1 
GOTO 2 
4 CONTINUE 
OPEN(UNIT =uout,FILE =’sa.out’»,STATUS =’unknown’) 
DO sa —a1e nl 
WRITE(uout, 105) 1,red(i),cred(i),blue(1),cblue(1) 
S CONTINUE 
n= n-1 
101 FORMAT(F10.7) 
102 FORMAT(A9) 
103 FORMAT(IS) 
104 FORMAT(6IS) 
105 FORMAT(I3,5(2x,F10.7)) 
RETURN 
END 
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SUBROUTINE evalfn(c,ecost) 


C** Sum of Square Differences (Linear Lanchester Model) 


201 


REAL _ cred(100),cblue(100),red(100),blue(100) 
INTEGER n 
COMMON /vars/cred,cblue,red,blue,n 


REAL ecost,r,b,c(4) 


ecost = 0.0 

DO 2011 =1,n 
r = (cred(i) - (c(1)*blue(i) + c(2)))**2 
b = (cblue(1) - (c(3)*red(i) + c(4)))**2 
ecost = ecost + r + b 

CONTINUE 

ecost = SORT (ecost) 

RETURN 

END 


67 


APPENDIX D 
GAMS REGRESSION MODEL 


$TITLE Least Squares Fit of Battle Data to Analytical Model 
$STITLE (Linear Lanchester Form - 2-stage fit) 


$OFFUPPER OFFSYMXREF OFFSYMLIST 
OPTIONS LIMCOL = 0, LIMROW = 0, SOLPRINT = OFF 
OPTIONS RESLIM = 50, ITERLIM = 10000, OPTCR = 0.0: 
faa coeaee DEFININIONS AND DA —————————— 


I increments aa a 


W weapons / 


RTANK 

RAPC 

BTANK 

BAPC 

BTOW /; 

PARAMETER FP(W) fire power / 

RTANK _ .920 
RAPC _ .780 
BTANK _ 1.000 
BAPC _ .850 


BTOW _ .750 /: 
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TABLE CHANGE(I,W) changes in wpn level in inc 


RTANK RAPC BTANK BAPC BTOW 


1 7 8 + 4 4 
2 8 8 3 4 + 
8 6 8 4 & + 
4 6 i 3, 3 ) 
» 6 6 Z , & 
6 6 6 3 3 6, 
i 5 6 Z Zz Z 
8 5 6 Z 2 3 
9 + 3 Us di i 
10 5 > 2 2 Z 
Ul 4 5 1 Z 1 
12 4 4 1 1 fies 


TABLE LEVEL(I,W) level of wpn at start of inc 


RTANK RAPC BTANK BAPC BTOW 
80 160 60 120 40 
(Sgts2en 26 6611660636 
65 144 33 112 32 
59 136 49 109 28 

129 46 106 25 

47 123 44 103 22 
Ahi aie 100 19 
BO oo 95 
31 105 37 £496 ~~ 14 

10 Za 1005.35 _ 94 aa 

Nal 222 95 335 OG all) 

WZ lon 90 32 190 on; 


Oo OIA NB WN 
SN 
QW 


PARAMETER RED(I) red firepower at increment 1; 


RED() = FP("RTANK") * LEVEL(,"RTANK") 
+ FP("RAPC") * LEVEL(,"RAPC"); 


PARAMETER CRED(I) change in red firepower in increment; 


CRED(I) = FP("RTANK") * CHANGE(I,"RTANK") 
+ FP("RAPC") * CHANGE(I,"RAPC’); 
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PARAMETER BLUE(I) blue firepower at increment 1; 


BLUE(I) = FP("BTANK") * LEVEL(I,"BTANK") 
aris FP("BAPC’) x LEVEL(I,"BAPC’) 
+ FP("BTOW") * LEVEL(I,"BTOW"); 


PARAMETER CBLUE(I) change in firepower of blue in increment; 
CBLUE(I)= FP("BTANK") * CHANGE(,"BTANK") 


+ FP("BAPC") * CHANGE(L"BAPC’) 
+ FP("BTOW") * CHANGE(L"BTOW’); 


*---------- MODEL---------- 
POSITIVE VARIABLES 
A coef on blue level in equation red 
C coef on red level in equation blue 
BL(I) residual in blue equation 
RD(1) residual in red equation; 
VARIABLES 
B additive noise term 
D additive noise term; 
VARIABLE 
COST sum of squared error; 
EQUATIONS 
OBJ minimize sum of squared error 
DRED 1(I) Sq error increment no noise term 
DBLUE1(I) Sq error increment no noise term 
DRED2(I) Sq error in increment with noise term 
DBLUE2(I) Sq error in increment with noise term; 
: >>>>>>>>>> MINIMIZE << <2 


OBJ.. COST =E= SUM(I, RD(I)) + SUM, BL()) ; 

: >>>>>>>>>> SUBJEC! 1G Ge 
DRED1(I).. RDU) =G= POWER((CRED(I) - ((A*BLUE(I)+8B))),2); 
DRED2(I).. RDU) =G= POWER((CRED(I) - ((A*BLUE(I)+B))),2); 


DBLUEI(I).. BL(I) =G= POWER((CBLUE(I) - ((C*RED(I) +D))),2); 
DBLUEQ2(I).. BL(I) =G= POWER((CBLUE(I) - ((C*RED(I) +D))),2): 
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MODEL LSQ1 /OBJ,DRED1,DBLUE1/; 
SOLVE LSQ1 USING NLP MINIMIZING COST; 


A.FX = A.L; 
ox = CL; 


MODEL LSQ2 /OBJ,DRED2,DBLUE2/; 
SOLVE LSQ2 USING NLP MINIMIZING COST; 


DISPLAY COST.L,A.L,B.L,C.L,D.L 


gk 


APPENDIX E 
PROBABILITY SURFACE COMPARISON CODE 
PROGRAM datmat 
INTEGER first(2),next(2, 100) 


REAL ~ ecdf(2,100) 
COMMON /list/first,next,ecdf 


INTEGER r1(20),r2(20), lower, higher 
REAL dmatrx(6,6,500),cdfdat(6,20),r(40) 
REAL plow, phigh, tdiff(6,6) 


CHARACTER* 12 cfile 
DOUBLE PRECISION dseed 


C dseed = 7561883.d0 
dseed = 829847.d0 
OPEN(UNIT=8,FILE =’cdfile.dat’, STATUS =’old’) 
n = 0 
111 CONTINUE 
n=n+1 
READ(8,211,END =S5S5) cfile 
OPEN(UNIT =9,FILE =cfile, STATUS =’old’) 
m = 0 
112 CONTINUE 
m=m+1 
READ(9,212,END =56) cdfdat(n,m) 
GOTO 112 
56 CONTINUE 
CLOSE(9) 
GOTO 111 
55 CONTINUE 
CLOSE(8) 
OPEN(UNIT= 10,FILE=’datmat.out’”»,STATUS =’unknown’) 
OPEN(UNIT = 11,FILE =’pdat.out’",STATUS =’unknown’) 


ve 


110 


Ly 


115 


116 


114 
1S 


119 
118 
117 
211 
Gs 
713 


DO 113i =1,5 
DO 114] = 1+1, 6 
lower = 0 
higher = 0 
PRINT™*, Comparing surfaces for : ’,1,’ vs ’,] 
DO 1101 = 1, 20 
ecdf(1,l1) = cdfdat(i,1) 
ecdf(2,1) = cdfdat(j,1) 
CONTINUE 
CALL tcdf(tdiff(i,j)) 
WRITE(11,*) 71,j’,1,j,’ tdiff. = ’,tdiff(i,j) 
PRINT *,’tdiff = ’,tdiff(i,j) 
DO 991 = 1, 20 
r(1) = cdfdat(i,l) 
r(1+20) = cdfdat(j,1) 
CONTINUE 
DO 1161 = 1,50 
CALL roll(20,40,r1,r2,dseed) 
DO 115 k = 1, 20 
ecdf(1,k) = r(r1(k)) 
ecdf(2,k) = r(r2(k)) 
CONTINUE 
CALL tcdf(dmatrx(1,j,1)) 
IF(dmatrx(,j,l) .LT. tdiff(i,j)) lower = lower +1 
IF(dmatrx(ij,l) .GT. tdiff(i,j)) higher= higher +1 
CONTINUE 
plow = FLOAT(lower)/50.0 
phigh = FLOAT(higher)/50.0 
CONTINUE 
CONTINUE 


DO 117i =1,5 
DO 118 j = i+1, 6 
DO 119k = 1, 50 
WRITE(10,213) dmatrx(i,j,k) 
CONTINUE 
CONTINUE 
CONTINUE 
FORMAT(A12) 
FORMAT(F10.7) 
FORMAT(F10.4) 


END 
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SUBROUTINE. tcdf(diff) 


INTEGER first(2),next(2,100) 
REAL _ ecdf(2,100) 
COMMON /list/first,next,ecdf 


REAL deltax,deltay,x,y,diff,ratio 


DO 33 ne = 1 ,2 
DO 32 n = 1, 10 
next(nc,n) = 0 
CALL ordcdf(nc,n) 
CONTINUE 
CONTINUE 
diff = 0.0 
deltax = 1.0 
deltay = 1.0 
DO 21 x = 50, 200, deltax 
Do 22 y = 50, 200, deltay 
ratio = (x/y)**2 
CALL cdf(1,ratio,prob1) 
CALL cdf(2,ratio,prob2) 
diff = diff + (abs(prob1-prob2) * (deltax*deltay)) 
CONTINUE 
CONTINUE 


FORMAT(F10.7) 
END 
SUBROUTINE ordcdf(nc,n) 
INTEGER count,test,ltest,ttest,maxcnt 
INTEGER first(2),next(2, 100) 
REAL ~ ecdf(2,100) 
COMMON /list/first,next,ecdf 


maxcnt = 100 
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IF (n .EQ. 1) THEN 
first(nc) = 
next(nc,1) = 
count = 1 
RETURN 

ELSEIF (count .EQ. maxcnt) THEN 
PRINT *, COVERFLOW ERROR >>> CDF LIST IS FULL << <’ 
STOP 

RIESE 
count = count + 1 
test = first(nc) 

201 CONTINUE 
IF (ecdf(ne,n) .LT. ecdf(nc,test)) THEN 

IF (test .EQ. first(nc)) THEN 
first(nc) = n 

ELSE 
next(ne,ltest) = n 

ENDIF 

next(nc,n) = test 

RETURN 

RIESE 

ttest = next(nc,test) 

IF (ttest EQ. 0) THEN 
next(nce,n) = 0 
next(nc,test) = n 
RETURN 

ELSE 
Itest = test 
test = ttest 
GOTO 201 

ENDIF 

ENDIF 
ENDIF 
RETURN 


END 


as 


400 
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SUBROUTINE cdf(nc,pt,prob) 


INTEGER first(2),next(2,100) 
REAL ~ ecdf(2,100) 
COMMON /list /first,next,ecdf 


INTEGER nc,test 
REAL _ pt,prob 


num = 0 
test = first(nc) 
CONTINUE 
IF (pt .LE. ecdf(nc,test)) THEN 
prob = FLOAT(num)/10.0 
RETURN 
EESE 
num = num + 1 
test = next(nc,test) 
IF (test EQ. 0) THEN 
prob = 1.0 
RETURN 
ENDIF 
ENDIF 
GOTO 400 
END 
SUBROUTINE roll(n,m,r,rc,dseed) 


INTEGER r(n),rc(n) 
DOUBLE PRECISION dseed 


DO? 1 Salm 
ri) = 0 
CONTINUE 
CALL unif(dseed,u) 
num = INT((u*m)+-.5) 
DO 1j=1,1 
IF (rg) .EQ. num) GOTO 3 
CONTINUE 
rai) = num 
CONTINUE 
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k = 1 
DO 6i=1,m 
DO4j=1,n 
IF (rq) .EQ. 1) GOTO 5 
CONTINUE 
re(k) = 1 
k=k+1 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE UNIF(DSEED,U) 
INTEGER I 

REAL U 

DOUBLE PRECISION DENOM,DSEED 
DATA DENOM/2147483647.D0/ 


DSEED = DMOD(16807.D0*DSEED,DENOM) 
U = DSEED/DENOM 


RETURN 
END 


T7 


LIST OF REFERENCES 


. Box, G.E.P., Hunter, W.G., and Hunter, J.S., Statistics for 
Experimenters: An introduction to Design, Data Analysis, and 
Model Building, pp. 93-106, John Wiley & Sons, Inc., 1978. 


. JANUS(A) Basic User’s Tutorial, TRAC-RDM-TD 90-1 (DRAFT), 
US Army TRADOC Analysis Command - Monterey, June 1990. 


. Ingber, L., Mathematical Companion of the Computer Models to 
Exercise Data: Comparison of JANUS(A) to National Training 
Center Data, Proposal to: US Army TRADOC Analysis Command, 
Fort Leavenworth, KS, 29 September 1988. 


. Shy, J., "Jomini," ed. Paret, P., Makers of Modern Strategy, 
pp. 143-185, Princeton University Press, 1986. 


. Taylor, J. G., Lanchester Models of Warfare, Volume I, 
p. 52, Operations Research Society of America, 1983. 


. Taylor, J. G., Lanchester Models of Warfare, Volume I, 
p. 53, Operations Research Society of America, 1983. 


. Sandefur, J. T., Discrete Dynamical Systems: Theory and 
Applications, pp. 213-222, Clarendon Press, 1990. 


. Hartman, J. K., "Lecture Notes in Aggregated Combat 
Modelling,", Naval Postgraduate School, 1985. 


. Bitters, D., "The Economics of the Battlefield," paper 
presented at the Naval Postgraduate School, Monterey, CA, 
November 1990. 


10. Barr, D., Weir, M. D., and Hoffman, J., "Evaluation of 


Combat," Naval Postgraduate School, Monterey, CA, 
unpublished. 


11. Metropolis, W., and others, "Equation of State Calculations 


by Fast Computing Machines,” Journal of Chemical Physics, 
v. 21, pp. 1087-1092, 1953. 


78 


—— 


12. 


:>: 


14. 


15. 


16. 


Kirkpatrick, S., Gelatt, C.D., and Vecchi, M. P., 
"Optimization by Simulated Annealing,” Science, v. 220, 
13 May 1983. 


Johnson, D. S., and others, "Optimization by Simulated 
Annealing: An Experimental Evaluation; Part I, Graph 
Partitioning,’ Operations Research, v. 37:6, p. 868, 
November-December 1989. 


Johnson, D. S., and others, "Optimization by Simulated 
Annealing: An Experimental Evaluation; Part II, The 
Traveling Salesman Problem," (in preparation). 


Brooks, D. G., and Verdini, W. A., "Computational Experience 
with Generalized Simulated Annealing Over Continuous 
Variables,” American Journal of Mathematical and Management 
Sciences, v. 8:3 & 4, 1988. 


Lewis, P. A. W., Orav, E. J., Simulation Methodology for 


Statistician, Operations Analysts, and Engineers, Volume I, 
Wadsworth & Brooks/Cole Advanced Books and Software, 1989. 


719 


BIBLIOGRAPHY 


Aarts, E., and Korst, J., Simulated Annealing and Boltzmann Machines, John Wiley & 
Sons, Inc., 1989. 


Bowman, M., Integration of the NTC Tactical Database and JANUS(T) Towards a Combat 
Decision Support System, Master’s Thesis, Naval Postgraduate School, Monterey, CA, 
1989. 


Brooke, A., Kendrick, D., and Meeraus, A., Gams: A Users Guide, The Scientific Press, 
1988. 


Giordano, F. R., and Weir, M. D., A First Course in Mathematical Modeling, Brooks/Cole 
Publishing Company, 1985. 


Golub, G. H., and Van Loan, C. F., Matrix Computations, 2d ed., The John Hopkins 
University Press, 1989. 


Ingber, L., "Very Fast Simulated Reannealing,", Mathematical Computer Modelling, v. 12, 
1989. 


Johnson, D. S., and others, "Optimization by Simulated Annealing: An Experimental 
Evaluation; Part II, Graph Coloring and Number Partitioning,” Operations Research, (to 


appear). 


Kirkpatrick, S., "Optimization by Simulated Annealing: Quantitative Studies," Journal of 
Statistical Physics, 1984. 


Naval Postgraduate School Report 55-79-029, Calculation of State Probabilities for a 
Stochastic Lanchester Combat Model, by L. Billard, November 1979. 


Taylor, J. G., Lanchester Models of Warfare, Volume IT, Operations Research Society of 
America, 1983. 


Tovey, C., "Simulated Simulated Annealing," , American Journal of Mathematical and 
Management Sciences, v. 8:3 & 4, 1988. 


Upton, S. C., Statistical Mechanics Model of Combat, Master’s Thesis, Naval Postgraduate 
School, Monterey, CA, March 1987. 


80 


. Professor Maurice D. Weir, Code MA/WC 


INITIAL DISTRIBUTION LIST 


No. Copies 


. Defense Technical Information Center 2 


Cameron Station 
Alexandria, VA 22304-6145 


. Library, Code 0142 Z 


Naval Postgraduate School 
Monterey, CA 93943-5002 


. Professor Donald Barr, Code MA/BA 2, 


Naval Postgraduate School 
Monterey, CA 93943-5004 


. Professor Patricia Jacobs, Code OR/JC 1 


Naval Postgraduate School 
Monterey, CA 93943-5004 


— 


Naval Postgraduate School 
Monterey, CA 93943-5004 


. Colonel Frank R. Giordano 1 


Department of Mathematical Sciences 
United States Military Academy 
West Point, NY 10996 


. Major James Hoffman 1 


TRAC-Monterey 
Monterey, CA 93943-0692 


. Captain John R. Green 3 


Department of Mathematical Sciences 
United States Military Academy 
West Point, NY 10996 


. Chairman, Mathematic Department, Code MA 1 


Naval Postgraduate School 
Monterey, CA 93943-5004 


81 


10. Chairman, Operations Research Department, Code OR 
Naval Postgraduate School 
Monterey, CA 93943-5004 


82 













































































































































at +" E% € + 
— i ee a Fm 
SS fine 4 004 held Of 24 eri EAs 
aie 5 rn rh of 0S Ff ut de 4 ot y re Soph Ast robe 2, 
4 i ST rel dicen LEAN Eile oh eat, eS oe ey re re bee aa 
ites, Pee eh a rr aah IP ot og. FO, 9. ote es © gog oH 64. 5 Lid ue ae 4: Biedods 4 2! tee 
“yd Val a er ei Seed gf hit pO ie a aS, f- [ohhe: a Mr) ahoh arto of bag Siutianrd 
ey. 7 4 Pe Be hs) 2 Le ee rd We Gals oF he Be ped we Bie , f 4, Gut f 
8 4d it ne, PP : pe BRE Met hg AL $l Bota. if ates “pil « jee Rint Pek ah ¢ Mb nite 
aithey ry y ee ee te of tgs. gegig pp ER AG AA gies eA MEY Wr are 
Snes oreo he ok ne Re ate NPof sof ENA tod oll taba gi hee aoe ts 
#, : je aA AT), ‘ops x CO Pt Bae 1 ol 0 a fy en Moet says dat £4 tH 
pore tyo hh pee ANE on Net MEY MOE Wins gocens Ape : ae wig ‘ sar ste gs ie “$ 
74) : Sot Aeil torah ret ee z uaa t oot 04 Hr aaah 
oe etree, SAOOK.OES iyi say, See © DATS OAD 38 oan a PL Ned Bt ERY Pa Frm 
SLO eres : Fe PAO a sre ea ar DOR Be Ng! Poss PEG 3 a Rie aN AY 
; is ‘Hae hae tLe, Gee A ers Seat of ar 
Manta b peel ret tate} 3 mHilttSmamign hy tp f bie waht bn me ere thefo ewe fred HU 
agin ao i prieyd an eee PP Pea Care i es f auhoia $4 PP eral DAS soda 
a DO. dprien ge OE th oh oF aa ‘ ng. cE NN BML ng np vet fapel BEE LASS tof okie oe ' 
MEM atesh Pine A ON Ao OAs on gene tA diab! WAI ge ao ge Police tds of en 
ERR Pe RECS pea ree LSE TORE RY Seay Be AI Lt aight cata 
AT ef ay areng ” = sae! iad TRS op Y EL aad ER Un Dei an Beas 48 
feta ai ey aaa etn a — dtr Gin henner fee Wiel BLg hie bok 
. et peg eO MAR We pgm tog 7. ‘i Bhar ol ae Led Racd ro La 
Lag RAED 8 Wisk fae) ABTS Utes welanae " ME hal 8 ot ts baked ba #4 j 
rire errs TR cre A. 063 8 2s, taets Po af fi P 4 atce LP 
Mg he He OF TAH estates 4 1p NTOS OB Besegieot sag eet eue ier tbl VS mms ie he tet ana gine 
Num h sen dive fe age: Pi Oh OH Ort san PAPA teh 4Abed, sotat Parner oO Bie hilly Beco he og 
Fal sahil met sions. meatal: gated of, PB LOA PAIS tad bi try, Mth Gis Ahabcl dss at SEA onDeb tke diod BRM, AEE oo J, 
nd atin Let Die diet afi NT RIM aati Od i chun, igs aft .2 “4 MR disted sic L oft ie #8 
RO rtenan: ee rae wap et Hhtt CARN AME ORE Sm She wegihot ca. Chee ihe tv iol Bam \%an 
WM ane ° FLOP Finn Boh v4, OFA P08 Pome gue # OE 37:0 = san $€sid aot os Mie bate ae A ‘ Fj 
Patna ee AR veho he maaie eee icine nner Ge Oe ottae ¢ ges bet A ated ha hat inecvenen, fed ed BP on 
LE MONO eh etna ot carne conan Meet Cobe Rs Ce So eee) Po hdl teak: eb, bagel £ 
it amnanogne, ot ep On Ba NS BER BE uk, ar iw M4, 4 alee, hotell 4.0 
Roe wennan, an peor vedi wre gn eee na APL mph sad pee te Ahead sas tegen LP ppbaelapey ‘ 
Le Piahemeree PRET pie ae) nr hanah een or So ee eee 9 SEF 492108. Wy oe op Hale tol ie, fy A he Mink, pdt Pag oat - 
enn HA 1 aos wisn lpontde Pent 18 Hs tcf wtsor wsyy Py PRA L Ls Et sid, eam : 
SEF Ab ach rivegee oes ge yom eK Mink Mth SEO, OF se BP OLS era AT Ripe $ aie = a amaecan oe i 
At NAM wr iteaesi ce weet Oe tegen g ahh lon oY For Peet rg s 
Pid Sete mes shee aM oh ang Fahd cc 0% Caner fo Se ee aa ; - Sa 
forcings ph nhhagheinen eT Oe habeehinta Com yl Satet fn rary mn e Othittce —— . 
mt an acarar” § Cana Sena rg Serre C21 aisha og rane bet nes Cah oS Pchel.or-rBahence fart ol ' 
aot eNe . : vr © St PAP hs Bde nlom a a ope liateea 6 rps mes" 
ne arene gute Pw pended 000Gb Andon 5 ww sec, nine: Pal horpenp Atpor eyelid! pipe P 
we arte ps ? Phat all? ba ren . fi 
ecg wre Pepa Pde Soe 2s lanteah a ack ve mr Y Pay ; ee FO ae og - 
Pees oe, Od neg Le OD sed 0 8 BF AF 06: tn) ee we Nip heeft F 
mpete a Bee ge PANN hah bad nee Cobh Lk ACU ide Sry. . 
ws SO PUN OE IN: Pgs » mledaa ne e t AP af < rs 
DIMA RB ee Oe Nate - i 
a ? 
s 
Mite 
te 
‘ir 
CF Tin, oo a 
Fan oA pe oeste Pe 
7 Gf o% * os * , 
A repr ad te ee 
‘ 8 ry a ‘ 
tev One Ma oe 
mM Dem ew. o& ¢ Fo 
wtih tos ta Mt 4 « # 
ida to oe F 
‘ ’ fam ¢ as f ute o.0 Fi 
ARO a Ae - erty *, : eh ates 
Metphe tA Peg bores + De spine 7 pr ae “ + 
CCR mE. 8 a a 7 fh ood tphen stn ae 
“eur * A OF aU ol Be 8 ay 88 ee IV el ok ta nag , 
ihe ds Leet Te Operetta eT e ane P ‘ s 
a AP eet Oy onal Ly a Rig on vase uw Py oe q 8 
A nue © ge 7 see " pgearse y- = 
OP Pore ¢ wana) Fes Or Oud 2 * ; : 4 
iad a led Ee "7" 9 - * 4 oe peed s 7 
THD oF ow. tin phahed cart ete see > J = ane 
* in@rror Wh [APY vt ee ap fs oo aia 
-. oo +e a 1° <4. mel om DP og ” "] - a “ee id 
“Fh wer ral me 8 ow 6 ath ett iT cela sz, oemea 9 oe pet ‘ 
Pr weed TBM yCCA. $50 a apes.. * pol # Po. os es , . 
Pt ep wey mms ¢ § Shanes ieie . ; ? ae 
Aare PVE cee 08 aren ; 3 : - Ae 
oe 8 matt Oeten, otat.a 4 @ ~ mn : p 
Awndaa. Phe on E ’ e ‘ 
s we ‘ Ae bus of . 
a - + of 
ogery "ie art ar oa ° 
i oe * a | € 
ow ¢.¢ ie ¢ ’ 
ow . ’ 
eee : e * ¢ 
ws ow = Poe as t 
4 : ’ U * 
= a ’ , s 
, fF pam ane a Pa Ps 
‘ ¢ feo , 
ayrer 4 St 4 ; x #5 ; 
inn Ge 
: fp os 
> # - : 
= Vee 
= y ° ae - — = 4 ms He * 
p ae i ye : ” ¢ ¢ , 
, i 4 
2 . . 
— . a tol 
= owt @i¢ o 
” , ] = ; 
. + : a 
= 
- rt) a 
rs 4 7 
= we + = ' 
: 
~ 6, 
hw 
hen 1G ~ 
4 me a oe —— 
’ ~~ =f & 
< r=" 3 
: weak, 
bd J 
] | - 2 * 
‘ e 
ae hl ba 
“ ~4% 3 
. 7 re = ay 2 4 
ea kp 
: ~ ae hae 
a 
‘ : 
~ Cif 
: " ‘ any 
/~ 
; =< rs * . 
; = e a 
* ‘ © 
= = pam ne me, & we S . *~ 
i SE Sens a8, 2H 4 
is “ — ~ . = = eee # oyee the 
a “= s 7. UMes ° 4 ye = > ‘ “ im, 
by ae, de be ee BO Secs ah ee 9 " 
< vo e * ae a ae) eo a, iy te, sabe mee BS =o ‘ 
a a deh whe, Fm ane te Men bye * 
Py 
P om “ 
~ e A e 
Wet « 
sath a's 
a oh 
~N 
— 
t 
% mh & 
" mas 4 
=. 7 ®& 
Ue nemeeni, ——— 6 x 
aos Saya *<% ’ , . 
Sipe yee ol to te he G0; ay, . 
we ‘ Seco F 
oye 6 g : 
abbot len be ba ee ~< 
SN Wherhas dln es aaa ¥ a i we 
rae 
4 
. os 
rae me 
Fe Seep ree | : 5 - ss wv 
NO Seth Fe Ny 0.4 Phe Eh ; : ae 
“PP RP 0G Sy oe by: : r 4 I ‘ 
: trie @ abe 
wi , Pager : . ie Res . 
ete a . ’ ? v ey = 4% ‘ 
inte tag ean Seg te aun 1. ea 
Se EA om Hyde 1 ta 979 *) tawe UV oF ‘ + 
See Se QA Eye SS Ly tre Meany MHS — Sirk, 
CN sales ing on lta ON fee hrf FS vee tale BR Nd Saareecs Pet pei 4a okt Mich dnind oe Ie : a 
WP UN eit, ‘ ye SATE eS ARAN QeK papi ie neck” Bie ESky, jew i ” 
Sban r hilerboitee telat tee Syst rt Made tn to St Le “an N tak ad . 
Mie eae Sh reindag) waa ete Um nbeiaeree sy, Sy BS saab bh ya) . ou te an | ot the LE Ng et? 5 ? ‘ 
Ree Te ee Ee ae Petal ven Sarge ie a TSS 
7 ‘ Rett ’ . ‘ on S 
Sy rN lita eee ns Re He Maps 9, cr, he te i Y “RA pete: SBA el eae . A 
heme red SAWS Dey my gOM Sarat hngtreennen MANN Swit t 
Ber eA a bistalene OLY Tetra bilan Cen tee = a Fe 6 
; ; by # PART Ce eye uw 
Soa eins the sey eet bite talent Yee Was , uhm e mms te hy So, “ ts 
uinderas : % A] Mery, | * Tey 
. x" wats ay, ; Py 
: Ab & ale Ge, ee ‘ ¥ 
ee a ome, wide 1S eee 
tire ee trecaegy ae So a ial ot Seat eee: re 
te Mish rep Ue Ake ie tee re 
nw RSE ULe ANE Gey Woveh-ane = Hp : UY 
MTR Ar ns arom te ce wee oe” VENK Oe . ‘ 
OW netieed Pe ¢ , 
Cor rts = 
to ths Se OVER YS Oey ? 
re rare J 
* = Peale Me by : e SNe ® = 
a bite et ed aa " Forng oi rey = 
ae SAE STN ws anpeden weeey ‘ rae tae” gt rr <a 
ity ies sae etm Stunna ebb rata een ry a B 
tae Le Ere IC eR Toe et iT é ™ : 
“ON @ Sane bpp D aT RRA 8 ORROPE OE 5. Une veh vee 
en’ arg gy Y° hod Mag, ie wy * 
Se eENe oaing Sere : ; 
Sen « 
‘ F , WV IN oy 2056 
ery, . ween baa Lak © 4 
ae! ate al beat eh erp ‘ 
: betenm th 1? 7 4 oe mice ss 
llamas te tk toe 



























« . + 
“bs & Rbhetei ae 4 . L a, +; . = 
PY. ep pe -aprcaay a5 KYSER ery ra r . 
BM TATE Ppp pee vy 997 8 wey y Sa des ue ke st 
Gees sek WPS hy Ys “> Wop TPPe ne @ a 4 v 
Ory WEP STU e 180 Sra wreer PAS ‘ 4 
Wd Ve PA yyy any “J hed AAD ai oe ee TEP ea a i 
o-Uv & eh ey VOY We? BIL + we we ¥ s ‘ 
Sy een ah enh St ae oiacy See 
<a . i ¥ 5 ah oy ee Sl Bowe =H ¥ VePey 
i ls woe WEE oy la fs > a bn ee 7 ~y 
icine cca creram athe chow arty satan ae ergs Sh 
ee WIL, Ve ns ‘ in FANG oR AS a 3 PNP EAN ly we > , 
WP ee thine Le TR irene wees PVG wanG pad. » Se sf ent 
x Se are ji A watem tre they att: Wo eae Cee ED Le Re p en 
Shon, ee et a8, o Yet ania, eatin: x SAR OR. Madd tnt SPR Ng ts Ly “fe, 
SAeoeeatasenres titan ek ee Nive, PS SP VP res opiy “arn by mane it are mh g's 
Pete eek. Cae 2 AN heroes Seta ta at ESAS ‘ 
Carico ay fos sare sreemyt tae Wigs rene. a Sala say aa gerd atoate om ne sols a aa a i i am 
y ) %?, «7 1h, he nt A sae 0 apdeatinn HS ae 4 
meee ahaceeertome ests te AF ee ARES OS O neeretiiteten outs 
NARA Bed eee hen, 





I 


i 


we 8 


my 





E 


Wied DUDL 

wl 
MH 

3 276 


| 
8 


KN 


| 
HAL 
00016428 9 


Hi 





‘ ‘ , ’ 
a ss . ’ 
, ‘ ‘ ’ 
v s 
# 
; ' 
t id 
’ 
2 Py ‘ 
i) 2 ‘ t . 
/ f 
4 U 4 . 
‘ ’ 
t on 
4 ‘ ‘ 
4 
‘ 
Y ‘ 
z ‘ 
‘ 
7 ‘ 
r 
. an) 
. 
‘ . : t 
. 
a 
. 
f 
6 
: 
‘ 
' 
. 
4 
° 
‘ ‘ 
‘ , 
° 
a 
‘ « ‘ 
’ 
‘ 
+ 
ry 
Pas 
” 
‘ 
é . 
t 
. 
. 
° 
: 
"ft 
' 
* 
‘ 
° 
‘s 
« 
* 
« 
x 
& 
so 
a 
d 
a 
. 
s 
r 
‘ 
* 
te 
F 
é 
‘i? 
3 © 


